function Opt = mysstateopt(This,Mode,varargin)
% mysstateopt  [Not a public function] Prepare steady-state solver options.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2015 IRIS Solutions Team.

% Run user-supplied steady-state solver:
% 'sstate=',@func
if length(varargin)==1 && isa(varargin{1},'function_handle')
    Opt = varargin{1};
    return
end

% Run user-supplied steady-state solver with extra arguments:
% 'sstate=',{@func,arg2,arg3,...}
if length(varargin)==1 && iscell(varargin{1}) ...
        && ~isempty(varargin{1}) ...
        && isa(varargin{1}{1},'function_handle')
    Opt = varargin{1};
    return
end

% Do not run steady-state solver:
% 'sstate=',false
if length(varargin)==1 && isequal(varargin{1},false)
    Opt = false;
    return
end

% Do run steady-state solve with default options:
% 'sstate=',true
if length(varargin)==1 && isequal(varargin{1},true)
    varargin(1) = [ ];
end

% `Mode` is either `'verbose'` (direct calls to `model/sstate`) or
% `'silent'`; the mode determines the default values for `'display='` and
% `'warning='`.
Opt = passvalopt(['model.mysstate',Mode],varargin{:});

%--------------------------------------------------------------------------

if This.IsLinear
    % Linear sstate solver
    %----------------------
    % No need to process any options for the linear sstate solver.
else
    % Non-linear sstate solver
    %--------------------------
    Opt = xxFixOpt(This, Opt);
    Opt = xxBlocks(This, Opt);
    Opt = xxDisplayOpt(This, Opt);
    Opt = xxOptimOpt(This, Opt);
    Opt = xxLogOpt(This, Opt);
    Opt = xxBounds(This, Opt);
end

end




function Opt = xxDisplayOpt(This,Opt) %#ok<INUSL>
if islogical(Opt.display)
    if Opt.display
        Opt.display = 'iter';
    else
        Opt.display = 'off';
    end
end
end




function opt = xxOptimOpt(this,opt) %#ok<INUSL>
% Use Levenberg-Marquardt because it can handle underdetermined systems.
oo = optimset( ...
    'display',opt.display, ...
    'maxiter',opt.maxiter, ...
    'maxfunevals',opt.maxfunevals,  ...
    'tolx',opt.tolx, ...
    'tolfun',opt.tolfun ...
    );
try %#ok<TRYNC>
    oo = optimset(oo, 'algorithm', {'levenberg-marquardt', 0.01});
end
if ~isempty(opt.optimset)
    if iscell(opt.optimset)
        opt.optimset(1:2:end) = regexprep(opt.optimset(1:2:end), '[^\w]', '');
        oo = optimset(oo, opt.optimset{:});
    else
        oo = optimset(oo, opt.optimset);
    end
end
opt.optimset = oo;
end




function opt = xxFixOpt(this, opt)
% Process the fix, fixallbut, fixlevel, fixlevelallbut, fixgrowth,
% and fixgrowthallbut options. All the user-supply information is
% combined into fixlevel and fixgrowth.
ixCanBeFixed = this.nametype==1 | this.nametype==2 | this.nametype==4;
list = { 'fix', 'fixlevel', 'fixgrowth' };
for i = 1 : length(list)
    fix = list{i};
    fixAllBut = [fix,'allbut'];
    
    % Convert charlist to cellstr.
    if ischar(opt.(fix)) ...
            && ~isempty(opt.(fix))
        opt.(fix) = regexp(opt.(fix),'\w+','match');
    end
    if ischar(opt.(fixAllBut)) ...
            && ~isempty(opt.(fixAllBut))
        opt.(fixAllBut) = ...
            regexp(opt.(fixAllBut),'\w+','match');
    end
    
    % Convert fixAllBut to fix.
    if ~isempty(opt.(fixAllBut))
        opt.(fix) = ...
            setdiff(this.name(ixCanBeFixed),opt.(fixAllBut));
    end
    
    if ~isempty(opt.(fix))
        fixPos = mylookupname(this,opt.(fix),[1,2,4]);
        ixValid = ~isnan(fixPos);
        if any(~ixValid)
            utils.error('model:mysstateopt', ...
                'Cannot fix this name: ''%s''.', ...
                opt.(fix){~ixValid});
        end
        opt.(fix) = fixPos;
    else
        opt.(fix) = [ ];
    end
end

fixL = false(1,length(this.name));
fixL(opt.fix) = true;
fixL(opt.fixlevel) = true;
fixG = false(1,length(this.name));
% Fix growth of endogenized parameters to zero.
fixG(this.nametype==43) = true;
if opt.growth
    fixG(opt.fix) = true;
    fixG(opt.fixgrowth) = true;
else
    fixG(:) = true;
end
% Fix optimal policy multipliers. The level and growth of
% multipliers will be set to zero in the main loop.
if opt.zeromultipliers
    fixL = fixL | this.multiplier;
    fixG = fixG | this.multiplier;
end
opt.PosFixL = find(fixL);
opt.PosFixG = find(fixG);
end




function opt = xxBlocks(this, opt)
% Run Blazer.
b = Blazer( );
b.Name = this.name;
b.Eqtn = this.eqtn;
b.IxOccur = this.occurS;
b.IxEndg = this.nametype==1 | this.nametype==2;
b.IxCanBeEndg = this.nametype==1 | this.nametype==2 | this.nametype==4;
b.IxEqtn = this.eqtntype==1 | this.eqtntype==2;
if ~isempty(opt.endogenise)
    b = endogenize(b, opt.endogenise);
end
if ~isempty(opt.exogenise)
    b = exogenize(b, opt.exogenise);
end
b = createBlocks(b, opt.blocks);
if b.IsSingular
    throw( exception.Steady.STRUCTURAL_SINGULARITY );
end
eqtnBlk = b.EqtnBlk;
nameBlkL = b.NameBlk;
nameBlkG = b.NameBlk;

% Finalize sstate equations.
eqtnS = myfinaleqtns(this, opt.growth);

nBlk = length(nameBlkL);
blkFunc = cell(1,nBlk);
ixAssign = false(1,nBlk);
% Remove variables fixed by the user.
% Prepare function handles to evaluate individual equation blocks.
for ii = 1 : nBlk
    % Exclude fixed levels, fixed growth rates, and refreshed values from
    % the list of optimised names.
    nameBlkL{ii} = setdiff(nameBlkL{ii}, opt.PosFixL);
    nameBlkL{ii} = setdiff(nameBlkL{ii}, this.Refresh);
    nameBlkG{ii} = setdiff(nameBlkG{ii}, opt.PosFixG);
    nameBlkG{ii} = setdiff(nameBlkG{ii}, this.Refresh);
    
    if isempty(nameBlkL{ii}) && isempty(nameBlkG{ii})
        continue
    end
    
    nameLPos = nameBlkL{ii};
    nameGPos = nameBlkG{ii};
    eqtnPos = eqtnBlk{ii};
    iiBlkEqtn = eqtnS(eqtnPos);
    
    % Check if this is a plain, single-equation assignment. If it is an
    % assignment, remove the LHS from `eqtn{1}`, and create a function
    % handle the same way as in other blocks.
    doTestAssign( );
        
    % Create a function handle used to evaluate each block of
    % equations or assignments.
    if true % ##### MOSW
        blkFunc{ii} = str2func(['@(x,dx) [',iiBlkEqtn{:},']']);
    else
        blkFunc{ii} = mosw.str2func(['@(x,dx) [',iiBlkEqtn{:},']']); %#ok<UNRCH>
    end
end

% Index of level and growth variables endogenous in sstate calculation.
ixEndgL = false(size(this.name));
ixEndgL([nameBlkL{:}]) = true;
ixEndgG = false(size(this.name));
ixEndgG([nameBlkG{:}]) = true;

% Index of level variables that will be always set to zero.
ixZeroL = false(size(this.name));
ixZeroL(this.nametype==3) = true;
if opt.growth
    ixZeroG = false(size(this.name));
    ixZeroG(this.nametype >= 3) = true;
else
    ixZeroG = true(size(this.name));
end

opt.EqtnBlk = eqtnBlk;
opt.NameBlkL = nameBlkL;
opt.NameBlkG = nameBlkG;
opt.IxAssign = ixAssign;
opt.BlkFunc = blkFunc;
opt.IxEndgL = ixEndgL;
opt.IxEndgG = ixEndgG;
opt.IxZeroL = ixZeroL;
opt.IxZeroG = ixZeroG;

return




    function doTestAssign( )
        % Test for plain assignment: One equation with one variable solved
        % for on the LHS.
        if length(iiBlkEqtn)>1 ...
                || length(nameLPos)>1 || length(nameGPos)>1
            return
        end
        namePos = nameLPos;
        if isempty(namePos)
            namePos = nameGPos;
        end
        lhs = sprintf('-(x(%g))',namePos);
        nLhs = length(lhs);
        % The variables that is this block solved for is the only thing on
        % the LHS but does not occur on the RHS.
        xn = sprintf('x(%g)',namePos);
        ixAssign(ii) = strncmp(iiBlkEqtn{1},lhs,nLhs) ...
            && isempty(strfind(iiBlkEqtn{1}(nLhs+1:end),xn));
        if ixAssign(ii)
            iiBlkEqtn{1}(1:nLhs) = '';
        end
    end
end




function opt = xxLogOpt(this, opt)
% xxLogOpt  Create the list of log-plus and log-minus levels,
% `Opt.IxLogPlus` and `Opt.IxLogMinus`, based on user options `'Unlog='`
% and `'LogMinus='`.

unlogList = opt.Unlog;
if ischar(unlogList)
    unlogList = regexp(unlogList,'\w+','match');
end
logMinusList = opt.LogMinus;
if ischar(logMinusList)
    logMinusList = regexp(logMinusList,'\w+','match');
end
conflict = intersect(unlogList,logMinusList);
if ~isempty(conflict)
    utils.error('model:mysstateopt', ...
        'This name is used in both ''Unlog='' and ''LogMinus='': ''%s''.', ...
        conflict{:});
end

% Positions of unlog variables.
unlogPos = mylookupname(this,unlogList,[1,2]);
ixValid = ~isnan(unlogPos);
if any(~ixValid)
    utils.error('model:mysstateopt', ...
        'This name cannot be used in ''Unlog='': ''%s''.', ...
        unlogList{~ixValid});
end

% Positions of log minus variables.
logMinusPos = mylookupname(this,logMinusList,[1,2]);
ixValid = ~isnan(logMinusPos);
if any(~ixValid)
    utils.error('model:mysstateopt', ...
        'This name cannot be used in ''LogMinus='': ''%s''.', ...
        unlogList{~ixValid});
end

% Create lists of log-plus and log-minus levels.
opt.IxLogPlus = this.IxLog;
opt.IxLogMinus = false(size(this.IxLog));
opt.IxLogMinus(logMinusPos) = true;
opt.IxLogPlus(logMinusPos) = false;
opt.IxLogPlus(unlogPos) = false;
opt.IxLogMinus(unlogPos) = false;
end




function opt = xxBounds(this, opt)
nName = length(this.name);
bndL = opt.levelbounds;
bndG = opt.growthbounds;
nBlk = length(opt.NameBlkL);
lowerL = cell(1,nBlk);
upperL = cell(1,nBlk);
lowerG = cell(1,nBlk);
upperG = cell(1,nBlk);
ixValidL = true(1,nName);
ixValidG = true(1,nName);
for iBlk = 1 : nBlk
    % Level bounds
    %--------------
    posL = opt.NameBlkL{iBlk};
    listL = this.name(posL);
    [lowerL{iBlk}, upperL{iBlk}, ixValidL] = ...
        doBlock(listL, posL, bndL, ixValidL);
    
    % Growth bounds
    %---------------
    posG = opt.NameBlkG{iBlk};
    listG = this.name(posG);
    lowerG{iBlk} = -inf(1,length(listG));
    [lowerG{iBlk}, upperG{iBlk}, ixValidG] = ...
        doBlock(listG, posG, bndG, ixValidG);
end

if any(~ixValidL)
    utils.error('model:mysstateopt', ...
        ['Lower and/or upper bound for steady-state level ', ...
        'of this variable has the wrong sign: ''%s''.'], ...
        this.name{~ixValidL});
end
if any(~ixValidG)
    utils.error('model:mysstateopt', ...
        ['Lower and/or upper bound for steady-state growth ', ...
        'of this variable has the wrong sign: ''%s''.'], ...
        this.name{~ixValidG});
end

opt.LowerL = lowerL;
opt.UpperL = upperL;
opt.LowerG = lowerG;
opt.UpperG = upperG;

return




    function [lowerB, upperB, ixValid] = doBlock(list, nameId, bnd, ixValid)
        nList = length(list);
        lowerB = -inf(1,nList);
        upperB = inf(1,nList);
        for jj = 1 : nList
            % Name of the variable.
            name = list{jj};
            % Position of the variables in This.name.
            isLogPlus = opt.IxLogPlus( nameId(jj) );
            isLogMinus = opt.IxLogMinus( nameId(jj) );
            try
                lo = bnd.(name)(1);
            catch
                if isLogPlus
                    lo = 0;
                else
                    lo = -Inf;
                end
            end
            try
                hi = bnd.(name)(2);
            catch
                if isLogMinus
                    hi = 0;
                else
                    hi = Inf;
                end
            end
            if isLogPlus
                lo = log(lo);
                hi = log(hi);
            elseif isLogMinus
                % Swap lower and upper bounds for log-minus variables.
                hi0 = log(-lo);
                lo = log(-hi);
                hi = hi0;
            end
            if imag(lo)==0 && imag(hi)==0
                % Assign only if both lower and upper bound is ok.
                lowerB(jj) = lo;
                upperB(jj) = hi;
            else
                % Report problem for this variables.
                ixValid( nameId(jj) ) = false;
            end
        end
    end
end
